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Abstract 


The  main  goal  of  the  research  project  is  to  use  sensitivity  analysis  to  design  actuators 
and  sensors  for  distributed  parameter  systems  (DPS).  The  choice  of  sensor/actuator  for 
a  control  problem,  its  action  and  placement  thereof,  is  addressed  by  computing  the 
sensitivity  of  the  state  variables  with  respect  to  small  changes  in  a  parameter  character¬ 
izing  placement.  Investigations  for  parameterized  sensor/actuator  placement  indicate 
that  computational  challenges  exist  for  certain  types  of  models.  Parameters  determining 
placement  tend  to  enter  the  PDE  models  in  ways  that  affect  the  differentiability  of  the 
state  variables.  There  can  be  a  loss  of  regularity  between  the  model  equations  and  that 
of  the  corresponding  sensitivity  equations.  When  using  Continuous  Sensitivity  Equation 
Methods  (CSEMs),  the  formulation  and  regularity  of  the  mathematical  model  governing 
the  underlying  system  behavior  must  be  well  understood.  Simple  models  involving  ac¬ 
tuator  placement  contain  discontinuous  coefficients,  and  the  explicit  parameter  location 
dependence  of  these  discontinuities  provides  an  a  priori  indication  that  sensitivity  vari¬ 
ables  lack  the  same  smoothness  properties  as  those  of  the  state.  The  issue  is  particularly 
important  for  accuracy  and  convergence  of  numerical  sensitivity  calculations. 
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1  Technical  Summary 


This  work  concerns  the  development  of  computational  tools  to  incorporate  into  the  design 
of  control  systems — control  law,  observers,  choice  and  placement  of  sensors  and  actuators — 
for  flexible  wing  MAVs.  Such  vehicles  have  been  the  focus  of  research  efforts  supported  by 
AFOSR  and  have  unique  capabilities  that  can  be  utilized  in  urban  missions.  Capabilities 
such  as  flight  control  via  wing  shaping  have  been  demonstrated,  and  it  is  known  that 
these  MAVs  can  fly  at  increased  angle  of  attack  without  stall.  However,  fundamental 
aerodynamics  of  such  vehicles  is  not  well-understood,  and  computational  studies  to  enhance 
the  of  control  systems  are  now  being  developed.  A  vital  contribution  to  that  effort  is  the 
choice  and  placement  of  sensors  and  actuators  used  within  the  control  regime. 

The  main  goal  of  the  research  project  is  to  use  sensitivity  analysis  to  design  actuators 
and  sensors  for  distributed  parameter  systems  (DPS).  The  choice  of  sensor/actuator  for  a 
control  problem,  its  action  and  placement  thereof,  is  addressed  by  computing  the  sensitivity 
of  the  state  variables  with  respect  to  small  changes  in  a  parameter  characterizing  placement. 
Preliminary  investigations  for  parameterized  sensor/actuator  placement  indicate  that  com¬ 
putational  challenges  exist  for  certain  types  of  models.  Parameters  determining  placement 
tend  to  enter  the  PDE  models  in  ways  that  affect  the  differentiability  of  the  state  variables. 
There  can  be  a  loss  of  regularity  between  the  model  equations  and  that  of  the  correspond¬ 
ing  sensitivity  equations.  When  using  Continuous  Sensitivity  Equation  Methods  (CSEMs), 
the  formulation  and  regularity  of  the  mathematical  model  governing  the  underlying  system 
behavior  must  be  well  understood.  Simple  models  involving  actuator  placement  contain 
discontinuous  coefficients,  and  the  explicit  parameter  location  dependence  of  these  discon¬ 
tinuities  provides  an  a  priori  indication  that  sensitivity  variables  lack  the  same  smoothness 
properties  as  those  of  the  state  variables.  The  issue  is  particularly  important  for  accuracy 
and  convergence  of  numerical  sensitivity  calculations. 

For  the  last  three  years,  my  collaborators  and  I  have  investigated  various  aspects  of 
the  use  of  sensitivity  analysis  for  locating  sensors  and  actuators.  Computational  models  are 
required  for  control  design  and  sensitivity  calculations  (for  sensor  and  actuator  selection  and 
placement).  The  focus  of  this  research  has  been  the  mathematical  analysis  and  development 
of  computational  tools  for  sensitivity  approximation  for  models  where  the  parameter  of 
interest  determines  sensor/actuator  placement.  Investigations  indicate  that  computational 
challenges  exist  for  certain  types  of  models,  and  one  important  indicator  of  such  challenges 
is  the  manner  in  which  the  parameters  appear  in  the  underlying  mathematical  models. 
Parameters  determining  placement  tend  to  enter  the  pde  models  in  ways  that  affect  the 
differentiability  of  the  state  variables  with  respect  to  small  changes  in  such  parameters.  At 
the  very  least,  there  may  be  a  loss  of  regularity  between  the  model  equations  and  that  of  the 
corresponding  sensitivity  equations.  In  Section  1.6,  we  discuss  a  proof-of-concept  example 
which  shows  how  sensitivity  computations  can  be  obtained  efficiently  and  can  be  validated 
if  the  scale  of  the  problem  allows.  In  Section  1.7,  an  Euler-Bernoulli  Beam  model  with 
patch  actuators  is  presented  in  order  to  illustrate  the  computational  challenges  that  may 
occur,  and  a  corresponding  steady  state  model  problem  is  also  examined  to  identify  some 
of  the  mathematical  issues. 

1.1  Methodology  and  Theory  Used  in  the  Research 

In  this  section,  an  overview  of  the  types  of  models  that  have  been  examined  in  the  research  is 
given.  Notation,  methodology  and  theory  is  also  outlined  and  will  be  referenced  throughout 
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the  report. 


1.2  Partial  Differential  Equation  Models  and  Their  Discretizations 

When  we  speak  of  a  “model” ,  we  are  referring  to  a  state  equation  that  describes  the  gov¬ 
erning  physics  of  the  problem  to  be  controlled.  Since  the  focus  of  this  report  examines 
actuator  placement,  we  omit  the  measurement  equation  that  models  the  sensor  dynamics. 
Specifically,  we  consider  the  simplified  underlying  model  considered  in  this  research  to  be  a 
linear  system  of  partial  differential  equations  of  the  form 

x(t)  =  Ax(t)  +  Bu(t),  x(0)  =  .To,  (1.1) 

where  x(t)  is  the  state  in  the  state  space  X  and  u{t)  is  the  control  effected  by  the  actuators. 
The  operator  A  represents  the  linear  dynamics  of  the  system,  and  B  represents  the  physics 
of  the  actuator  input  to  the  system.  Depending  on  the  application,  the  state  vector  may 
be  composed  solely  of  flow  velocities,  or  perhaps  velocities  and  structural  components  if  we 
are  considering  the  coupling  between  the  vehicle  and  flow.  In  a  more  general  discussion, 
one  would  include  several  important  terms  in  the  model,  most  importantly,  a  nonlinear 
operator  describing  the  nonlinear  dynamics  of  the  system  to  be  modelled.  Other  terms 
such  as  a  disturbance  term  that  can  be  used  to  incorporate  noise  or  unmodeled  dynamics 
as  well  as  a  measurement  term  that  provides  a  measurement  of  the  state  provided  by  the 
sensors  should  also  appear.  However,  for  the  discussion  given  in  this  report  we  simplify 
the  governing  equation  to  a  linear  model  and  focus  on  the  most  interesting  aspects  of  the 
sensitivity  analysis  research  that  has  been  conducted  for  this  project. 

This  type  of  model  will  be  used  as  guidance  in  algorithm  development.  For  example, 
theory  for  sensitivity  analysis  for  distributed  parameter  systems  uses  the  PDE  model  as  the 
starting  point  [17,  69,  71,  72,  74].  Similarly,  design  of  convergent  control  designs  relies  on 
properties  of  the  underlying  PDE  system. 

Deriving  a  finite  dimensional  model,  one  can  obtain  a  discretized  version  of  (1.1)  of  the 
form 

xN(t)  =  ANxN{t)  +  BNuN{t ),  xN(0)  =  x$  (1.2) 

where  N  is  related  to  the  grid  used  in  discretization,  and  the  state  vector  xN(t )  is  obtained 
by  applying  an  appropriate  approximation  scheme  such  as  a  finite  element  discretization. 


1.3  System  Sensitivities 

In  this  research  project,  we  make  extensive  use  of  the  sensitivity  of  the  state  of  the  system, 
or  control,  to  a  parameter  of  interest,  a.  The  parameter  may  represent  a  sensor  or  actuator 
location,  its  length,  or  some  other  characteristic  that  would  affect  the  overall  closed  loop 
system  performance  and  or  robustness.  If  we  denote  the  sensitivity  of  the  state  variable 


xQ(t) 


dx(t) 

da 


then  we  can  derive  the  sensitivity  equation  for  the  system  defined  in  (1.1)  by  implicitly 
differentiating  to  yield 


xa(t)  =  Axa(t)  +  Aax(t )  +  Bua(t)  +  Bau(t), 

xQ(0)  =  0. 


(1.3) 
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Issues  with  respect  to  existence,  regularity,  and  calculation  of  the  sensitivities  for  partic¬ 
ular  systems  provide  fundamental  questions  for  research.  Controller  sensitivities  will  be 
presented  after  the  overview  of  the  basic  control  method  given  in  the  following  section. 

1.3.1  The  Linear  Quadratic  Regulator  Problem 

Given  a  system  defined  by  any  of  the  models  described  by  equations  (1.1)  or  by  (1.2),  the 
associated  linear  quadratic  regulator  (LQR)  problem  is  as  follows: 

Determine  the  control  u(t)  in  an  admissible  set,  U,  that  minimizes  the  cost  functional 

/»oo 

J{u)=  /  (Qx,  x)x  +  (Ru,  u)u  dt, 

Jo 

subject  to  the  constraint  that  x(t)  in  state  space  X  satisfies  the  state  equation,  and  Q  and 
R  are  state  and  control  weighting  operators,  respectively.  Assuming  full  state  feedback, 
optimal  control  theory  gives  the  feedback  control  in  the  form 

u(t)  =  -RTlB*Ylx(t)  =  —Kx(t),  (1.4) 

where  K  is  the  gain  operator,  and  II  is  the  solution  of  the  control  algebraic  Riccati  equation 

A*U  +  UA-UBR~1B*U  +  Q  =  0.  (1.5) 


Functional  Gains 

For  many  PDE  control  problems,  the  control  law  (1.4)  can  be  written  as 

u(t)  = —[Kx(t)\  =  /  k(s)x(t,  s)ds,  (1.6) 

Jo. 

for  example  see  [19,  20,  48,  50].  In  the  case  that  there  is  more  than  one  actuator  input,  e.g., 
more  at  more  than  one  spatial  location,  £,  there  will  be  a  functional  gain  for  each  location 
and  the  inputs  can  be  represented  as 

ut(t)  = -[Kx(t) ](£)  =  f  k^(s)wc(t,s)ds.  (1.7) 

Jo 

The  kernels  k^(s)  are  called  functional  gains,  and  have  been  utilized  for  reduced  order 
controller  design  and  sensor  placement  in  [4,  5,  21,  41,  49,  53]. 


1.4  Control  System  Sensitivities 

Computing  sensitivities  of  the  control  law  with  respect  to  parameters  is  one  of  the  primary 
research  interests.  If  we  define  the  control  sensitivity 


ua(t) 


du 

da 


then  the  associated  sensitivity  system  is  given  by  equation  (1.3)  and 


ua(t)  =  —R~1  [B*Hxa(t)  +  (R*na  +  R*n)x(t)] , 


where  na  solves  the  Lyapunov  equation 

A*Ua  +  UaA  =  Q ,  (1.8) 
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with 


A  =  A  -  BR~lB*Yl  and  Q  =  -A*aU  -  n Aa  +  n (BR~lB*a  +  BaR~lB*)B  -  Qa. 

This  Lyapunov  equation  is  derived  by  implicitly  differentiating  through  the  Riccati  equation 
(1.5)  with  respect  to  a.  Note  the  use  of  the  notation  Aa  =  for  the  sensitivity  of  the  A 
operator,  and  similarly  for  Ba  and  Qa.  We  note  that  in  general,  the  control  weighting,  R, 
does  not  depend  on  a. 


1.5  Results  of  Funding 

1.6  Investigations  on  Sensitivity  Computations  for  the  Heat  Equation 
with  Respect  to  Actuator  Placement 

This  section  addresses  an  actuator  design  problem  for  the  simple  ID  heat  equation.  The 
idea  is  that  the  choice  and  placement  of  sensors  and  actuators  for  a  control  problem  could 
be  addressed  by  computing  the  sensitivity  of  the  state  of  the  system  with  respect  to  a 
characteristic  of  the  actuator  placement.  To  illustrate  the  ideas,  we  consider  the  sensitivity 
of  the  state  variable,  temperature  in  this  case,  to  actuator  location.  This  work  provides  the 
preliminary  computations  that  demonstrate  proof-of-concept  for  the  combination  of  CSEMs 
with  basic  optimal  control  techniques.  The  sensitivity  of  the  temperature  with  respect  to 
actuator  placement  is  investigated  using  a  Continuous  Sensitivity  Equation  Method,  see 
[17,  43,  71,  72,  74,  75]  for  examples.  Consider  the  heat  equation  model  for  the  situation 
when  an  actuator  location  is  specified  by  the  parameter  a: 


wt(t,  £)  =  ew#(t,  0  +  b(£;  a)u(t) 


(1.9) 


w(t,  0)  =  0,  w(t,  1)  =  0 

u;(0,£)  =  w0(£)  (1.10) 

with 

6(£;  a)  = 

Defining  the  state  variable  x(t)  =  w(t,  •)  and  using  standard  PDE  theory,  one  can  formulate 
(1.9)-(1.10)  in  the  form  of  (1.1)  posed  in  an  appropriate  Sobolev  space.  We  forego  the 
presentation  of  these  details  for  the  sake  of  brevity.  However,  we  note  that  the  B  operator 
above  depends  explicitly  on  the  parameter  of  interest  a;  this  fact  is  important  for  the 
sensitivity  analysis  that  follows. 

Many  of  the  computational  issues  related  to  sensitivity  analysis  depend  heavily  on  the 
correct  mathematical  formation  for  the  problem.  Here  we  outline  the  application  of  the 
Continuous  Sensitivity  Equation  Method  (CSEM)  to  the  LQR  control  problem  for  the  heat 
equation  given  above.  The  following  discussion  uses  the  temperature  (state)  sensitivity,  the 
n  operator  sensitivity  and  the  controller  sensitivity,  respectively,  defined  as 


wa(t,Q 


dw(t,£) 

da 


nn 


dn 

da  ’ 


ua{t) 


du(t ) 
da 


To  derive  the  appropriate  sensitivity  equations,  one  implicitly  differentiates  equations  (1.9)- 
(1.10)  with  respect  to  the  parameter  a.  The  resulting  sensitivity  equation  with  boundary 
and  initial  conditions  is  given  by 


(■ wa)t(t ,  £)  =  e(wa)^(t,  £)  +  2(£  -  a)b(£;  a)u(t)  +  6(£;  a)uQ(t)  (1.11) 
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wa(t,  0)  =  0,  wa(t,  1)  =  0 

Wa{  0,0  =  0. 

Further,  we  note  that  this  equation  can  be  posed  in  a  simplified  form  of  equation  (1.3)  as 

xa(t)  =  Axa(t)  +  Bua(t)  +  Bau(t) ,  xQ(0)  =  0  (1.12) 

The  state  feedback  control  system  is  given  by  equations  (1.1),  and  the  optimal  feedback 
control  by  (1.4)  determined  by  the  algebraic  Riccati  equation  (1.5).  The  control  sensitivity 
can  be  computed  as  described  in  Section  1.4.  In  this  example,  both  operators,  B  and 
Ba,  have  integral  representations,  and  it  is  an  easy  calculation  to  see  that  the  B  operator 
varies  in  a  differentiable  manner  with  respect  to  the  actuator  location  parameter  a.  In  the 
following  section  we  outline  the  discretization  and  sensitivity  approximations  for  both  the 
state  variable  as  well  as  the  control. 


1.6.1  Numerical  Experiments 


Here,  a  brief  summary  of  the  numerical  experiments  for  the  actuator  placement  modelis 
given.  The  initial  condition  used  for  the  numerical  calculations  is  the  piecewise  constant 
function 


™o(£)  =  w(0,£) 


-1,  £<=[0,0.5] 
l,  £e  (0.5, 1.0], 


and  the  thermal  diffusivity  coefficient  for  all  of  the  numerical  results  presented  here  is  given 
by  e  =  0.01.  The  state  equation  in  (1.12)  is  discretized  using  a  standard  finite  element 
method  in  the  spatial  coordinate  with  cubic  B-spline  basis  functions,  and  a  time  solver 
from  Matlab’s  suite  of  ODE  solvers  is  used  for  the  time  stepping.  As  this  is  a  fairly 
standard  treatment  of  the  discretization  and  simulation,  we  omit  details.  Matlab’s  lqr 
and  lyap  functions  are  used  to  solve  the  algebraic  Riccati  equation  to  approximate  n  and 
to  approximate  the  Riccati  sensitivity,  na,  respectively. 

In  Figure  1,  an  example  of  a  temperature  approximation  using  piecewise  cubic  basis 
functions  is  shown  for  the  case  where  a  =  0.7  and  N  =  40  basis  functions  are  used.  Figure 
2  contains  the  corresponding  temperature  sensitivity  for  those  same  parameter  values.  Note 
that  the  sensitivity  gives  the  qualitative  information  that  one  would  expect  for  the  case  when 
a  =  0.7.  That  is,  the  behavior  of  wa  in  Figure  2  shows  that  the  temperature  w  is  most 
sensitive  in  the  region  around  £  =  0.7,  and  that  the  temperature  is  most  sensitive  for  early 
values  of  t  before  the  diffusion  and  the  controller  have  driven  the  temperature  to  be  very 
near  equilibrium.  The  associated  control  and  control  sensitivity  computations  are  given  in 
Figure  3  and  4.  The  behavior  of  these  functions  is  also  as  expected.  The  control  effort 
is  large  initially  and  quickly  decays  to  0.  Similarly,  the  control  is  most  sensitive  to  small 
changes  in  a  for  small  values  oft.  As  the  t- variable  increases,  the  control  is  less  sensitive  to 
changes  in  a.  In  the  following  section,  we  include  a  brief  discussion  on  the  effort  to  validate 
the  sensitivity  calculations  for  this  simple  example. 


1.6.2  Validation  of  Sensitivity  Calculations 

This  section  briefly  demonstrates  that  for  the  simple  example  of  the  ID  heat  equation, 
one  can  validate  the  sensitivity  calculations  shown  in  the  previous  section.  The  follow¬ 
ing  shows  that  this  can  be  done  for  the  temperature,  the  control  variable,  or  the  Riccati 
operator  approximations.  To  clarify,  the  numerical  approximations  of  the  sensitivities  for 
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Figure  1:  Temperature  Approximation  w(t ,  £)  with  N  =  40  and  a  =  0.7 


various  variables  are  compared  with  finite  difference  calculations.  For  example,  the  control 
sensitivity,  ua  is  compared  with  a  centered  finite  difference  approximation  so  that 


Ua  (* ; «) 


uN(t;  a  +  Aa)  —  uN (t;  a  —  A  a) 
2Aa 


The  superscript  notation  uN  refers  to  the  fact  that  the  calculation  is  based  on  a  finite  el¬ 
ement  discretization  in  space,  and  for  the  results  presented  below,  a  reasonably  fine  mesh 
of  N  =  40  is  used  in  order  to  guarantee  an  accurate  approximation  for  the  variables  of 
interest.  Hence,  the  convergence  that  we  present  here  refers  to  a  convergence  of  the  finite 
difference  approximation  as  one  allows  A  a  — >  0.  Table  1  gives  results  for  each  of  the  sensi¬ 
tivity  variables  wa,  uQ,  and  Ua  for  a  specific  value  of  the  a  parameter.  Similar  results  were 
observed  for  a  range  of  parameter  values.  Although  sensitivity  calculations  obtained  via 
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0.005 


Figure  3:  Control  Approximation  u(t)  with  N  =  40  and  a  =  0.7 


Figure  4:  Control  Sensitivity  Approximation  ua(t)  with  N  =  40  and  a  =  0.7 


finite  difference  techniques  can  be  difficult  to  obtain  for  large  scale  problems,  it  is  appro¬ 
priate  to  validate  our  sensitivity  calculations  whenever  possible.  For  the  case  of  this  simple 
model  problem,  sensitivity  calculations  obtained  by  the  CSEM  approach  are  validated  by 
comparison  to  centered  finite  difference  sensitivity  computations.  In  the  following  section, 
one  will  see  why  such  a  validation  is  not  always  feasible  for  certain  types  of  sensitivity  vari¬ 
ables,  particularly  for  certain  types  of  sensitivity  variables  related  to  actuator  location  for 
some  mathematical  models.  An  investigation  of  actuator  location  for  a  boundary  control 
problem  with  the  2D  heat  equation  can  by  found  in  [32]. 


9 


Step  Size 

State  (cr  =  0.7) 

Control  (a  =  0.7) 

Riccati  (a  =  0.2) 

Act  =  10"1 
A  a  =  10-2 
Acr  =  10~3 
Acr  =  10-4 
Acr  =  10-5 

2.0600  x  10"2 
2.0782  x  10-5 
1.0914  x  10“6 
1.0900  x  10-7 
1.0900  x  10“8 

3.1782  x  10"4 
3.1802  x  10-7 
3.1801  x  10”10 
3.1418  x  10-13 
1.4816  x  10~14 

3.1000  x  10“3 
3.0808  x  10-5 
3.0808  x  10“7 
2.9993  x  10-9 

Table  1:  Error  Calculations  Showing  Convergence  of  Finite  Difference  Calculations  for  wa, 
ua  and  na. 

1.7  Computational  Issues  with  Sensitivity  Calculations  for  Actuator  Place¬ 
ment  on  an  Euler-Bernoulli  Beam 

In  this  section  we  discuss  an  actuator  placement  problem  for  a  cantilevered  Euler-Bernoulli 
beam  with  piezoceranric  patch  actuators  as  in  [8,  9,  10].  We  consider  the  sensitivity  of  the 
state — position  and  velocity  of  the  beam — to  patch  location  and  to  patch  length.  This  work 
provides  the  preliminary  computations  which  would  then  be  used  for  optimal  selection  of 
patch  length  and  location. 

The  framework  that  is  currently  being  developed  combines  CSEMs,  optimization  and 
LQR  control  techniques  to  gain  insight  into  sensor  and  actuator  placement  for  control  design 
for  physical  systems  governed  by  a  partial  differential  equations.  However,  this  section  shows 
that  the  sensitivity  calculations  for  actuator  placement  can  present  many  mathematically 
interesting  and  numerically  challenging  problems. 

This  discussion  relates  preliminary  sensitivity  calculations  for  a  cantilevered  Euler- 
Bernoulli  beam  model  with  a  pair  of  piezoceranric  patch  actuators  placed  on  either  side 
of  the  beam  at  the  same  specified  spatial  location  (see  Figure  5).  We  consider  the  case 
when  the  patches  are  excited  out-of-phase  which  results  in  pure  bending  of  the  beam.  As 
in  the  previous  example,  the  sensitivity  of  the  displacement  of  the  beam  with  respect  to 
patch  placement  is  investigated  using  the  CSEM.  We  incorporate  both  Kelvin- Voigt  and 


Figure  5:  Diagram  of  beam  with  a  pair  of  piezoceranric  patches  of  length  L  located  on  either 
side  of  the  beam  on  the  interval  x  £  [a,  a  +  L\. 


linear  viscous  (air)  damping.  The  length  of  the  beam  is  £,  the  height  is  given  by  b  and  the 
thickness  is  denoted  as  h.  The  length  of  the  patch  is  given  by  L,  and  the  patch  location 
is  given  by  the  interval  [a,  a  +  L\.  Let  w(t,x)  denote  the  displacement  of  the  beam  at 
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time  t  and  position  x.  The  motion  of  the  cantilevered  beam  is  governed  by  the  PDE  with 
boundary  conditions 


pAwtt  T  divwt  T  (dfcvlwtxx'jxx  T  (EIwxx)xx  —  gi^i  x)>  (1.13) 

w(t,  0)  =  w^(t,  0)  =  0  (1.14) 

EIwxx(t,  tj  +  dkvIwtxx(t,  £)  =  0,  (1.15) 

(EIwxx)x(t,t)  +  (dkvIwtxx)x(t,£)  =  0.  (1.16) 

and  with  the  initial  deflection  and  initial  velocity  denoted  by 

w(0,  x)  =  wq(x),  wt( 0,x)=wi(x).  (1.17) 


The  coefficients  all  represent  material  properties  of  the  beam  at  a  certain  spatial  location  x: 
p  is  the  mass  density,  A  is  the  cross-sectional  area,  I  is  the  moment  of  inertia,  E  is  Young’s 
modulus,  and  div  and  dj~v  are  the  coefficients  of  air  and  Kelvin- Voigt  damping,  respectively. 

The  presence  of  the  patches  results  in  discontinuities  in  these  coefficients,  which  can  be 
expressed  as  follows 


pA(x) 

=  pAi  +  pA2[Hi(x)  -  H2(x)], 

(1.18) 

dkyd (x) 

=  dkvh  +  dkvh[Hi(x)  -  H2(x)}, 

(1.19) 

EI(x ) 

=  Eh  +  EL2[Hi(x)  -  H2(x)}. 

(1.20) 

The  notation  H\{x)  denotes  the  Heaviside  function  taking  on  the  value  of  0  for  x  £  [0,  a] 
with  jump  discontinuity  located  at  x  =  a,  the  left  end  of  the  patch,  and  H2  denotes  the 
heaviside  function  with  jump  discontinuity  located  at  x  =  a  +  L,  the  right  end  of  the 
patch.  The  constants  pA\,EI\ ,  and  dkvh  correspond  to  the  density,  flexural  rigidity,  and 
Kelvin- Voigt  damping  properties  of  the  beam,  while  the  constants  pA2,EI2 ,  and  dkvh 
correspond  to  those  of  the  patch,  respectively.  The  patches  influence  the  beam  by  exerting 
a  moment  force  on  that  section.  Therefore  the  spatial  influence  of  the  control  is  described 
by  a  difference  of  Heaviside  functions,  and  the  control  function  for  the  beam  is  given  by 


g(t,x)  =  k[Hi(x)  -  H2(x)]xxu(t). 


(1.21) 


The  constant  n  is  a  parameter  characterizing  the  patch  properties,  and  u(t)  is  the  voltage 
applied  to  the  patch  at  time  t.  For  a  more  thorough  treatment  of  the  model  development 
and  the  specific  forms  of  the  beam  and  patch  parameters,  see  [10]. 

One  should  observe  that  in  the  beam  equation  (1.13),  the  damping  ( dkvIwtxx)xx ,  stiff¬ 
ness  ( EIwxx)xx ,  and  control  term  g(t,x)  all  contain  spatial  derivatives  of  the  Heaviside 
functions;  moreover,  the  actuator  location  parameter  is  embedded  in  the  definition  of  these 
Heaviside  functions.  Consequently,  the  PDE  should  be  mathematically  interpreted  using 
the  variational  (or  weak)  formulation.  The  weak  form  is  also  convenient  for  the  numerical 
simulations  which  use  the  finite  element  method  for  approximation  in  space.  The  variational 
formulation  of  the  beam  equation  is  given  in  [10] 

r£  rl 

/  pA(x;a)wtt{t,x)(j)xx(x)  dx  +  /  dkvI{x\  a)wtxx(t,  x)<f>xx(x)  dx 
Jo  Jo 

r£  r£ 

-I-  /  EI(x;  a)wxx(t,  x)cf>xx(x)  dx  =  /  g(t,  x\  a)<j)xx{x)  dx,  (1-22) 
Jo  Jo 
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where 


g(t,x;a)  =  k[Hi(x)  —  H2(x)]u(t). 

We  use  this  formulation  below  to  discuss  computational  difficulties  with  numerically  ap¬ 
proximating  the  sensitivities. 


1.7.1  Computational  Difficulties  with  the  Sensitivity  Equation 

The  parameter  of  interest  for  this  study  is  the  parameter  determining  the  patch  location, 
a,  and  we  consider  both  the  displacement  and  the  control  term  to  depend  explicitly  on  a. 
We  denote  this  dependence  as  w(t,  x)  =  w(t,  x\  a).  The  sensitivity  of  the  state  with  respect 
to  the  patch  location  is  defined  by 

.  .  dw  ,  . 

wait,  x;  a)  =  ——w(t,x;a). 
oa 


In  [73] ,  the  PI  numerically  approximated  wa  by  deriving  a  variational  sensitivity  equation. 
This  equation  was  obtained  by  differentiating  through  the  above  variational  equation  (1.22) 
for  the  state  with  respect  to  a.  This  procedure  causes  derivatives  of  the  Heaviside  and  delta 
functions  to  appear  in  the  sensitivity  equation.  The  regularity  of  the  sensitivity  equation 
can  play  an  important  role  in  determining  the  appropriate  computational  treatment  of  the 
problem. 

To  illustrate  the  main  difficulties  with  the  beam  sensitivity  equation,  we  remove  the 
time  dependence  to  focus  on  the  issue  of  regularity  in  the  spatial  domain.  We  also  remove 
one  spatial  discontinuity  from  the  piecewise  constant  functions  to  further  simplify  the  pre¬ 
sentation.  These  modifications  allow  us  to  study  a  similar  problem  which  has  a  closed  form 
solution  for  the  state  and  sensitivity.  Therefore,  instead  of  the  variational  equation  (1.22), 
consider  the  following  steady  problem:  find  w  G  V  =  {u  G  H2(0,£):  u(0)  =  0,  ux(0)  =  0} 
such  that 

r£ 

/  EI(x;a)wxx(x)(t>xx(x)dx  =  /  g(x]  a)4>xx{x)  dx  (1-23) 

Jo  Jo 

for  all  (j)  £  V.  The  functions  EI{x ;  a)  and  g(x\  a)  are  piecewise  constant  and  are  given  by 


EI(x ;  a) 


Eli,  0  <  x  <  a 
EI2,  a  <  x  <  i 


g(x-,a) 


gi,  0  <  x  <  a 
<72,  a  <  x  <  £ 


This  problem  can  be  posed  in  an  abstract  differential  form;  however,  like  the  time  dependent 
beam  problem,  it  is  best  understood  in  weak  form  due  to  the  discontinuities  in  the  coefficient 
functions  El  and  g. 

We  assume  Eli,  EI2  >  0  which  implies  the  problem  has  a  unique  solution  by  the  Lax- 
Milgram  Theorem.  It  can  be  checked  that  the  solution  is  given  by 

,  \  _  /  bix2/2,  0  <  x  <  a  .  . 

w^x,a  |  b2X2 /2  +  a{bi  —  b-2)x  —  a2{bi  —  b2) /2,  a  <  x  <  £  ’ 


where  6*  =  gi/  El  j  for  i  =  1,2.  The  solution  is  differentiable  with  respect  to  a  and  the 
sensitivity  wa(x;a)  is  given  by 


f0,  0  <  x  <  a 

\  (61  —  b2)x  —  a(bi  —  62),  a  <  x  <  £ 


(1.25) 
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The  sensitivity  is  continuous  and  differentiable  in  x  with  derivative 


(wa)x(x;  a) 


0,  0  <  x  <  a 

61  —  62,  a  <  x  <  i 


If  61  /  62,  then  (wa)x  is  not  differentiable  and  therefore  wa  ^  V.  Therefore,  the  sensitivity 
does  not  have  the  same  spatial  regularity  as  the  state. 

This  observation  has  implications  for  computing  the  sensitivities.  We  can  proceed  in 
the  same  manner  as  described  above  for  the  beam  problem  and  formally  differentiate  the 
steady  problem  (1.23)  with  respect  to  a  to  obtain  the  sensitivity  equation 

EI(x-  a)(wa)xx(x]  a)4>xx(x)  dx  +  (Eh  -  EI2)iuxx(a)(fxx(a)  =  (gi  -  g2)(f)xx(a).  (1-26) 

Recall  that  1 j>  is  in  H2(0,£ )  and  so  the  pointwise  evaluations  <fxx(ot)  do  not  make  sense  in 
general.  Also,  a  direct  computation  shows 


wxx(x) 


6i/2,  0  <  x  <  a 
62/2,  a  <  x  <  l 


so  that  the  term  wxx(a )  appearing  in  the  sensitivity  equation  does  not  have  meaning  when 
61  h  b2.  In  particular,  the  function  w  possesses  a  second  derivative  only  in  the  weak  sense 
and  not  in  the  classical  sense  of  differentiability. 

Despite  these  observations,  we  can  follow  the  approach  taken  in  [73]  and  numerically 
approximate  the  sensitivity  using  the  formally  derived  sensitivity  equation  (1.26).  We  use 
the  finite  element  method  for  the  spatial  discretization  with  cubic  B-spline  basis  functions. 
Since  these  basis  functions  have  continuous  second  derivatives,  the  second  derivative  evalu¬ 
ations  4>xx(a)  are  well  defined.  To  approximate  wxx(a)  (which  does  not  exist  in  the  classical 
sense),  we  use  the  second  derivative  of  the  finite  element  solution  of  the  steady  problem. 
The  numerical  solutions  of  the  formal  sensitivity  equation  for  32  equally  spaced  nodes  is 
compared  with  the  exact  sensitivity  in  Figure  6. 

The  numerical  solution  of  the  formal  sensitivity  equation  is  not  close  to  the  true  sensi¬ 
tivity.  Furthermore,  as  the  mesh  is  refined  this  error  does  not  decrease  and  the  numerical 
sensitivity  does  not  converge.  This  may  not  be  surprising  since  this  approach  uses  cubic 
polynomial  basis  functions  to  approximate  a  function  whose  second  derivative  is  not  smooth. 
Hence,  the  accuracy  of  the  sensitivity  computations  and  the  reliability  of  the  Euler-Bernoulli 
beam  system  simulations  is  brought  into  question  by  the  results  of  this  steady  problem.  The 
discontinuous  coefficients  in  equations  (1.22)  and  (1.23)  along  with  the  explicit  parameter 
dependence  of  these  discontinuities  is  the  first  indication  that  such  numerical  issues  must  be 
confronted  for  this  problem.  This  steady  problem  shows  that  approximating  a  sensitivity 
using  an  ill-posed  sensitivity  equation  may  lead  to  large  errors. 


Remark:  It  is  interesting  to  note  that  there  is  a  direct  connection  between  the  formal 
sensitivity  equation  and  automatic  differentiation  (AD).  In  our  computations,  we  used  the 
same  number  of  finite  element  basis  functions  for  the  state  equation  and  the  formal  sensitiv¬ 
ity  equation.  In  this  case,  the  discretized  sensitivity  equation  can  be  derived  by  implicitily 
differentiating  the  discretized  state  equation  with  respect  to  the  parameter  a.  This  is  ex¬ 
actly  how  one  uses  AD  to  numerically  approximate  the  sensitivity.  Therefore,  the  widely 
used  AD  procedure  also  fails  to  give  accurate  sensitivity  approximations  for  this  problem. 

In  future  work,  the  PI  plans  to  investigate  how  to  correctly  formulate  the  sensitivity 
equation  in  order  to  obtain  numerically  reliable  sensitivities.  See  Section  ??  for  more  de¬ 
tails. 
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Figure  6:  Finite  element  solution  using  cubic  B-splines  (with  32  equally  spaced  nodes)  of 
the  formal  sensitivity  equation  (1.26)  compared  with  the  exact  sensitivity  (1.25). 


1.8  Conclusions  and  Future  Work 

The  heat  equation  example  in  Section  1.6  shows  that  sensitivity  information  has  potential 
to  provide  valuable  information  for  optimizing  sensor  and  actuator  location.  However,  the 
steady  beam  equation  example  in  Section  1.7.1  shows  very  clearly  that  when  using  the  con¬ 
tinuous  sensitivity  equation  method  to  compute  sensitivities,  the  formulation  and  regularity 
of  the  mathematical  model  governing  the  underlying  system  behavior  must  be  well  under¬ 
stood.  One  goal  of  this  research  effort  is  to  advance  mathematical  and  computational  tools 
to  accurately  compute  parameter  sensitivities  for  systems  governed  by  partial  differential 
equations. 

For  the  steady  equation  in  (1.23),  and  similarly  for  the  beam  example  in  (1 . 13)-(1. 17) , 
the  discontinuous  coefficients  as  well  as  the  explicit  parameter  dependence  of  these  dis¬ 
continuities  provide  the  researcher  with  an  upfront  indication  that  there  may  be  a  loss 
of  regularity  when  moving  from  the  original  state  equation  to  the  derivation  of  the  cor¬ 
responding  sensitivity  equation.  This  is  noted  specifically  in  equation  (1.26).  Although 
this  phenomenon  occurs  in  a  simple  model  problem,  this  issue  is  expected  to  appear  in 
more  realistic  problems,  such  as  the  placement  of  actuators  on  a  flexible  wing.  When  sen¬ 
sor/actuator  placement  is  parameterized  in  such  a  manner,  the  formulation  and  analysis  of 
the  corresponding  sensitivity  equations  must  be  done  carefully  to  insure  that  any  regularity 
issues  are  accounted  for.  As  discussed  in  Section  1.7.1,  the  regularity  issue  is  of  particular 
importance  for  accuracy  and  convergence  of  numerical  sensitivity  approximations. 

For  the  beam  problem  with  the  discontinuous  coefficients,  research  indicates  that  there 
are  many  existing  mathematical  options  that  can  be  used  to  help  address  the  formulation 
of  the  sensitivity  equation  as  well  as  the  numerical  computations.  In  particular,  one  can  use 
a  very  weak  formulation  of  the  sensitivity  equation  to  approximate  the  sensitivity  variables 
that  may  lack  the  same  smoothness  as  their  original  state  variable  counterparts.  This 
approach  has  been  successfully  implemented  in  [23]  to  approximate  sensitivities  with  respect 
to  interface  location  in  a  heat  conduction  problem.  Using  this  technique  avoids  the  regularity 
issues  encountered  in  the  formally  defined  sensitivity  equation  presented  in  equation  (1.26). 
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Another  potential  method  to  treat  the  regularity  issue  of  the  sensitivity  equations  in¬ 
volves  smoothing  the  discontinuous  coefficients  in  the  original  mathematical  model  prior 
to  the  sensitivity  equation  derivation  and  applying  any  discretization  technique.  However, 
preliminary  research  using  such  a  technique  indicates  that  there  are  computational  issues 
involved  in  accurately  approximating  the  state  gradients,  and  the  smoothness  of  the  sensi¬ 
tivity  is  still  in  question.  This  is  the  subject  of  a  paper  which  is  currently  in  preparation. 

In  addition,  one  can  investigate  the  use  of  the  immersed  finite  element  method  [57] 
in  order  to  obtain  numerically  reliable  sensitivities.  Other  approaches  such  as  discontin¬ 
uous  Galerkin  finite  elements  and  Nitche’s  method  are  also  possible.  In  any  case,  when 
sensor/actuator  location  is  in  question,  one  has  to  treat  the  derivation  of  the  CSEM  and 
corresponding  sensitivity  calculations  using  techniques  which  account  for  loss  of  smoothness 
in  such  sensitivity  variables. 

One  other  note  is  that  if  the  sensitivity  calculations  were  to  be  carried  out  using  the 
very  basic  finite  difference  techniques,  such  as  those  used  in  Section  3.1.2,  the  lack  of 
smoothness  of  the  sensitivity  variable  would  impact  the  accuracy  and  convergence  rates 
of  those  computations.  In  addition,  other  sensitivity  techniques  such  as  the  popular  and 
well-advanced  Automatic  Differentiation  (AD)  software  package  also  suffers  from  similar 
problems.  Regardless  of  which  of  the  popular  techniques  one  uses  for  the  sensitivity  analysis, 
the  mathematics  of  the  sensor /actuator  location  problems  must  be  handled  carefully. 
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Stokes  Equations,”  J.  R.  Singler,  submitted  to  Journal  of  Mathematical  Analysis  and 
Applications. 
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3  Technology  Assists,  Transitions  or  Transfers 


None 

4  Accomplishments  and  Successes 

The  research  has  shown  that  sensitivity  information  has  the  potential  to  provide  valuable 
information  for  optimizing  sensor  and  actuator  location  for  control  design.  However,  investi¬ 
gations  into  certain  types  of  model  problems  also  has  shown  that  when  using  the  continuous 
sensitivity  equation  method  to  compute  sensitivities,  the  formulation  and  regularity  of  the 
mathematical  model  governing  the  underlying  system  behavior  must  be  well  understood. 
Various  mathematical  models  have  been  investigated  with  varying  degrees  of  success.  Struc¬ 
tural  models  such  as  the  outlined  in  the  report  include  the  Euler-Bernoulli  beam  model 
with  discontinuous  coefficients  representing  actuator  placement  along  with  some  models 
that  arose  due  to  simplification  of  this  model.  Dr.  Faranak  Pahlevani  investigated  various 
aspects  sensitivity  analysis  with  respect  to  the  eddy  viscosity  parameter  in  Eddy  Viscos¬ 
ity  Models  which  are  a  subclass  of  Large-Eddy  Simulation  models  for  computational  fluid 
dynamics  applications.  Dr.  John  Singler  has  developed  a  theoretical  framework  for  the 
derivation  of  continuous  sensitivity  equations  for  a  certain  class  of  parabolic  PDE  models. 
In  addition,  he  has  collaborated  with  the  PI  on  various  aspects  of  sensitivity  computation 
for  the  Euler-Bernoulli  beam  model.  Dr.  Singler  has  also  investigated  the  use  sensitivity 
analysis  to  explain  and  predict  transition  to  turbulence  for  Burgers’  Equation  models  with 
the  goal  of  future  application  to  the  study  of  transition  to  turbulence  for  Navier-Stokes 
equations. 

The  successes  of  this  research  effort  include  the  advance  of  the  mathematical  and  com¬ 
putational  tools  for  accurate  computation  of  parameter  sensitivities  for  systems  governed 
by  elliptic,  parabolic  as  welll  as  nonlinear  elliptic  partial  differential  equations.  Through  the 
use  of  such  available  methods  as  the  Petrov-Galerkin  Finite  Element  Method,  accurate  sen¬ 
sitivity  computation  for  elliptic  interface  problems  were  shown  to  be  a  viable  option.  Other 
types  of  methods  including  Discontiuous  Galerkin  as  well  as  Immersed  Interface  methods 
are  also  very  likely  to  prove  useful  for  future  efforts  in  this  area.  Furthermore,  the  research 
has  clearly  shown  that  relying  on  black-box  solvers  for  accurate,  as  well  as  numerically 
convergent,  sensitivity  calculations  may  lead  to  problems.  The  research  has  clearly  shown 
that  for  mathematical  models  where  the  parameter  of  interest  corresponds  to  locations  of 
sensors  or  actuators,  the  derivation  and  simulation  of  sensitivity  equations  must  be  handled 
with  caution.  The  main  successes  of  the  project  include  the  discovery  of  examples  for  which 
standard  continuous  sensitivity  equation  methods  work  well  and  may  be  done  somewhat 
automatically  as  well  as  the  identification  of  certain  types  of  mathematical  models  where 
the  casual  sensitivity  equation  derivation  will  likely  lead  to  misleading  sensitivity  informa¬ 
tion.  This  may  ultimately  lead  to  inaccurate  sensitivity  information  used  for  the  placement 
of  sensors  and  actuators  in  control  systems  which  may  ultimately  lead  to  poor  performance 
of  various  control  systems.  Further  investigations  into  several  of  these  ideas  is  warranted, 
and  the  PI  and  collaborators  continue  to  investigate  many  of  the  ideas  uncovered  during 
the  life  of  the  funding. 
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